Current hysteresis and memory effect in a molecular quantum dot with strong 

electron-vibron interaction 



(N 
O 

o 

(N 

o 

Q 

oo 



13 

o 
o 



> 

(N 

O 

-i— > 



i 

c 

o 
o 



X 



A.S. Alexandrov 1,2 and A.M. Bratkovsky 1 
1 Hewlett-Packard Laboratories, 1501 Page Mill Road, 1L, Palo Alto, California 94304 
2 Department of Physics, Loughborough University, Loughborough LE11 3TU, United Kingdom 

(December 16, 2002) 

Theory of current hysteresis for tunneling through a molecular quantum dot (MQD) with strong 
electron-vibron interactions and attractive electron-electron correlations is developed. The dot is 
modeled as a d-fold degenerate energy level weakly coupled to the leads. The effective attractive 
interaction between polarons in the dot results in a "switching" phenomenon in the current- voltage 
characteristics when d > 2, in agreement with the results for the phenomenological negative-f/ 
model. The degenerate MQD with strong electron-vibron coupling has two stable current states in 
certain interval of the bias voltage below some critical temperature. 
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I. INTRODUCTION 

Strongly correlated transport through mesoscopic sys- 
tems with repulsive electron-electron interactions has re- 
ceived considerable interest in the past (see, for exam- 
ple, |l],H,H]J]|,|l,0,U ) ' an d continues to be the focus of 
intense experimental and theoretical investigation |^|,|l0| . 
The Coulomb interaction suppresses tunneling for certain 
range of applied voltages, leading to what is commonly 
called the Coulomb blockade. There is now a growing 
interest in molecular nanowires and quantum dots used 
as the "transmission lines" [ pT]JT^ ] and active molecular 
elements 1 13 14|] in molecular-scale electronics Jl3[ . 

A few experimental studies fL4| provide evidence for a 
molecular switching effect, when the current-voltage (I- 
V) characteristics show two branches with high and low 
current for the same voltage. The effect exists in simple 
molecules too |l5|]. It is important to identify the actual 
mechanism of switching. 

Recently we have proposed a negative— U Hubbard 
model of a d-fold degenerate quantum dot, with an intrin- 
sic non-retarded switching mechanism when d > 2 [ pH) . 
We argued that the attractive electron correlations could 
be caused by a strong electron-phonon (vibron) interac- 
tion in the molecule, and/or by the valence fluctuations. 

It has been recently demonstrated that the low-bias 
conductance of molecules is dominated by resonant tun- 
neling through coupled electronic and vibration levels 
|jr7f . Different aspects of the electron-phonon/ vibron (e- 
ph) interaction effect on the tunneling through molecules 
and quantum dots (QD) have been studied by several au- 
thors ||Jl||^f|,||,||,||Q. In particular, Glazman 
and Shekhter, and later Wingreen et al. |Q presented the 
exact resonant-tunneling transmission probability fully 
taking into account the e-ph interaction on a nondc- 
generate resonant site. Phonons produced transmission 
sidebands but did not affect the integral transmission 



probability. Li, Chen and Zhou |19|] studied the con- 
ductance of a double degenerate (due to spin) quantum 
dot with Coulomb repulsion and the e-ph interaction. 
Their numerical results also showed the sideband peaks 
and the main peak related to the Coulomb repulsion, 
which was decreased by the e-ph interaction. Kang pcj ] 
studied the boson (vibron) assisted transport through a 
double-degenerate QD coupled to two superconducting 
leads and found multiple peaks in the I-V curves, which 
originated from the singular BCS density of states and 
the phonon sidebands. Ermakov |2l]] calculated the I-V 
curves of a four-fold degenerate dot including both the 
onsite Coulomb and e-ph interactions. He obtained a 
switching effect in the numerical I-V curves, similar to 
that in the negative-?/ Hubbard model discussed by us 
recently p6| . However, using a Hamiltonian averaged 
over the phonons, Ermakov missed all phonon sidebands 
and obtained an unphysical population (n = 1) of each 
QD state. More recently Gogolin and Komnik |2f| have 
analyzed a nondegenerate QD, d = 1, coupled with a sin- 
gle phonon mode. They found a switching effect in the 
Born-Oppenhcimer approximation similar to that in our 
negative-?/ Hubbard model, but surprisingly for a non- 
degenerate case with d = 1. However, we have to mention 
that the Born-Oppenhcimer approximation does not ap- 
ply to a non-degenerate level, since there are no "fast" 
(compared to the characteristic phonon time 1/loq) elec- 
tron transitions within the dot. Despite differences in 
the models and approximations, Refs. E^ , plL^5| pointed 
to a novel mechanism of the volatile molecular memory 
caused by the e-ph or any other attractive electron cor- 
relations. 

Here we develop the analytical theory of a correlated 
transport through a degenerate molecule quantum dot 
(MQD) fully taking into account both Coulomb and e-ph 
interactions. We show that the phonon sidebands signif- 
icantly modify the switching behavior of the I-V curves 
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in comparison with the negative-?/ Hubbard model [ p.6| . 
Nevertheless, the switching effect is robust. It shows up 
when the effective interaction of polarons is attractive 
and the state of the dot is multiply degenerate, d > 2. 



II. STEADY CURRENT THROUGH MQD 

We apply the Landauer-type expression for the steady 
current through a region of interacting electrons, derived 
by Meir and Wingreen as (in units h = ks = 1) 

/oo 
du\fi(u>)-h(u)]TmTr I»G«H 
-oo L 



where /i( 2 )M 



exp[(^ + A (+) eV/2)/T] + 1 



(1) 



T 



is the temperature, A is the position of the lowest unoc- 
cupied molecular level with respect to the chemical po- 
tential. r(w) depends on the density of states (DOS) 
in the leads and on the hopping integrals connecting 
one-particle states in the left (1) and the right (2) leads 
with the states in MQD, Fig. 1. This formula includes, 
by means of the Fourier transform of the full molecular 
retarded Green's function (GF), G^(u>), the e-ph and 
Coulomb interactions inside MQD and coupling to the 
leads. Since the leads are metallic, electron-electron and 
e-ph interactions in the leads, and interactions of elec- 
trons in the leads with electrons and phonons in MQD 
can be neglected. We are interested in the tunneling 
near the conventional threshold, eV = 2 A, Fig.l, within 
a voltage range about an effective attractive potential \ U\ 
caused by phonons/ vibrons (see below). 

The attractive energy is the difference of two large in- 
teractions, the Coulomb repulsion and the phonon medi- 
ated attraction, of the order of leV each. Hence, \U\ is 
of the order of a few tens of one eV. We neglect the en- 
ergy dependence of T(u>) ~ T on this scale, and assume 
that the coupling to the leads is weak, T <C \U\. In this 
case G^fjjj) does not depend on the leads. Moreover 
we assume that there is a complete set of one-particle 
molecular states \fi), where G^(u>) is diagonal. With 
these assumptions we can reduce Eq.(j^) to 



I{V)=h dw[/i(«)-/ 2 (w)]pH, 



(2) 



allowing for a transparent analysis of essential physics 
of the switching phenomenon. Here Iq — eT and the 
molecular DOS, p(w), is given by 
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FIG. 1. Schematic of the energy levels and phonon 
side-bands for molecular quantum dot under bias voltage V 
(eV/2A = 0.75) with the coupling constant j 2 = 11/13. The 
level is assumed to be 4-fold degenerate (d = 4) with ener- 
gies A + rU, r = 0, . . . , (d — 1) (thick bars). Thin bars show 
the vibron side-bands with the size of the bar proportional 
to the weight of the particular contribution in the density 
of states (see text) in the case of one vibron with frequency 
ojo/A = 0.2 at T = 0. Only the bands in the energy win- 
dow (eV/2, —eV/2) (shown) will contribute to current at zero 
temperature. 



.(*) = 



e iHt c„e- lHt , 



9(t) = 1 for t > and zero oth- 



erwise. We calculate p(u>) exactly in Section III in the 
framework of the Hamiltonian, which includes both the 
Coulomb U c and e-ph interactions as 



H = e n +-V U°,i 



+ 53 K w i^m d i + H - c -) + Yl u M d i + V2). (4) 

Here e fi are one-particle molecular energy levels, — 
cJjCp the occupation number operators, c m and d q anni- 
hilates electrons and phonons, respectively, u q are the 
phonon (vibron) frequencies, and ^ m are e-ph coupling 
constants (q enumerates the vibron modes). This Hamil- 
tonian conserves the occupation numbers of molecular 
states n . Hence it is compatible with Eq. (0) . 



(3) 



where G^(lu) is the Fourier transform of G^(t) — 
—i0(t) ({^(i), ct}) , {••■,••■} is the anticommutator, 
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III. MQD DENSITY OF STATES 



We apply the canonical polaron unitary transformation 



e 1261, integrating phonons out, as 



H = e s He~ s , 



where 



S = ~ r V ilM d q ~ H - c ] 



(5) 



(6) 



is such that = —S. The electron and phonon operators 
are transformed as 



and 



dg = d q-J2 



(7) 

(8) 



respectively. Here 



Xn = exp 



^qdq - Ha 



The Lang-Firsov canonical transformation shifts ions to 
new equilibrium positions with no effect on the phonon 
frequencies. The diagonalization is exact: 



* = £< 



where 



Y^u q {d\d q + 1/2) + - Y EWW> 



(9) 



(10) 



where now electron and phonon operators are averaged 
over the quantum state of the transformed Hamilto- 
nian H. There is no coupling between polarons and vi- 
brons in the transformed Hamiltonian, and the electron 
and phonon averages are independent. The Heisenberg 
phonon operators evolve as 



(13) 



so that we find after thermodynamic averaging of the 
phonon correlator over phonon occupation numbers (us- 
ing the Weyl's identity for exponential operators) 



exp 



£ 



I 7m I 



q Sinh 2 



cos \ tot + i 



PUq 



cosh 



(3uj q 



, (14) 



where (3 — 1/T. Repeating the calculations for 
(X^X^t)) we find a simple useful relation 



(xlx^t)) = (x,{t)xlY 



(15) 



At low temperatures T <C uj q the phonon correlator ((Tj) 
simplifies to 



(X„(t)Xl) = exp 



£l7 M rV^-l) 



(16) 



Next, we introduce the multi-particle GFs, which will 
necessarily appear in the equations of motion for 

(c p (t)cl), as 

G<r>+)(t) = -i6(t) Wn ft A 



and 



is the interaction of polarons comprising their interac- 
tion via molecular deformations (vibrons) and non-vibron 
(e.g. Coulomb repulsion) U^,. To simplify the discus- 
sion, we shall assume, without losing generality, that the 
Coulomb integrals do not depend on the orbital index, 
i.e. U^> = U. 

The molecular energy levels are shifted by the polaron 
level-shift due to a deformation well created by polaron, 



(11) 



Applying the same transformation in the retarded GF we 
obtain 



G p (t) = -*0(t) <{ C/1 (t)* p (t), ctxt}) 

= -ie(t)[(c^t)cl)(x^t)xl) 

+ ( C +e M (t))(X+X M (i))], 



(12) 



1=1 



Then, using the equation of motion for the Heisenberg 
polaron operator, 



dt 



" £ 



(19) 



we derive the following equations for the multi-particle 
GFs, 



dG ( ;< +) (t) 



dt 



S(t)(l-n,) yi n 



A'1t^M2#---A i »=1 

+ [£ +(r- 1)U]G^+Ht) + UG% +1 >+\t), (20) 



and 



3 



dG 



(*) 



r-1 



(If 



= 5{t)n, II 



Tin 



+% + (r - l)U]G^-\t) + UG^-\t), (21) 

where = (c^c^ is the expectation number of electrons 
on the molecular level \x. 

We can readily solve this set of coupled equations for 
MQD with one (Mold degenerate energy level and with 
the e-ph coupling 7 Mg = 7 9 , which does not break the 
degeneracy. At zero bias voltage the empty level will lie 
by some energy A above the Fermi levels of the electrodes 
(Fig. 1). Assuming that = n, Fourier transformation 
of the set yields 



G( 1 '+)H = (l-n) 



Z r (n) 



d-l 

^ u - rU + i(5 ' 

r=0 



d-l 



where 6 = +0, and 



Z r (n) 



r=0 



ui — rU + iS 



(22) 



(23) 



(24) 



A. MQD Green's function at low temperatures 

(T < LU q ) 

It is the easiest to find the total Green's function at 
low temperatures. Indeed, transforming back to real time 
and using Eqs.^), @ and @, @ we arrive at 



d-i 



irUt 



G„(t) = -i6(t)Z ^Z r {n) 

r=0 

x (1 - n) exp ( l7w 
L V q 

+nexp l^|7 M | 2 e^* 



(25) 



where Zq = exp — |7 g | 2 J . This is an exact solution 
with respect to correlations and e-ph interactions which 
satisfies all sum rules. Expanding the exponents in the 
brackets we find its Fourier component as 



d-i 



G^(uj) = Z o y^Z r (n) 

■r=0 

oo _^ 

+ Eyi X hn-Jii 



i 



lv — rU + id 



1=1 qi,—qi 



u-rU + Yl k=1 < 



1 -n 



iS 



(26) 



If the e-ph interaction is weak, \ j q \ <§; 1, the essential con- 
tribution comes only from the first (phonon-less) term, 
and we recover the result of the Hubbard model [fill 



Z r (n) 



lu — rU + i5 



(27) 



At finite \j q \ > 1, the phonon side-bands become im- 
portant in Eq.(p6]), which is obviously in the form of the 
multi-polaron expansion. If one neglects the correlations, 
U = 0, a standard polaron GF p7| is recovered: 



= Z Q 



w + iS ^11 

1=1 qi,...qi 



1 — n 



J2k=i 



T, k =i UJ v k +i8 



(28) 



where we have applied the sum rule 

d-l 



E^r(n) = l. 



(29) 



r=Q 



B. MQD Green's function at finite temperatures, 
single vibron mode 

By applying the same method, as in the case of T = 0, 
and going over back to real time with the use of Eqs.(p2|), 
(^H) and (|l^), (|l5| ) we arrive at 



d-l 



-irUt 



G M (t) = -i9(t)ZJ2 z r(n)e 

r=0 

x (l- «) exp fx; l7w 

( expfx^ 



,, sinh^ 



COS Ldt + i — - 



COS Ojt — I 



fJUJg 



where now 



exp 



X |7g| 2 COth 



2 



(30) 



(31) 



In approximation, where we retain a coupling to a sin- 
gle mode with the characteristic frequency ujq and 7 9 = 7, 
we can expand the exponents in the temporal Green's 
function ( |30| ) in powers of exp (cut + i^^j ■ It is then 

trivial to find the Green's function in the frequency do- 
main as 



d-l 



G^)=Z^2z r (n)J2l l (0 



1=0 



4 



e 2 



+(l-6 l0 )e 2 ( 



lu — rU — Iujq + iS lu — rU + Iujq + iS 
1 — n 



lu — rU + lu>o + id 



lu — rU — Iujq + id 



(32) 



where £ = \j\ 2 / sinh ^p, 7; (£) is the modified Bessel 
function, and <5jfc is the Kroneker symbol. At low tem- 
peratures, where (3uiq > 1, ( < 1 and Ii (£) w //!, 
this expression gives (Eq) in the form 



ItI 1 



i=0 



1 — n 



Z! Vu; — r?7 — Ilu + ^ 



cj — rU + /c^o + ?'<5 



(33) 



The molecular DOS is readily found as an imaginary 
part of Eq.@: 



d-l 



r=0 



2=0 



e 2 [(1 — u)5(lu — rU - Iuq) + nS(uj - rU + Iluq)} 



+ (1 - Si )e 2 [nS(uj ~ rU - Ilu ) 
+(1 - n)5(u) - rU + lu Q )] 



(34) 



The important feature of the DOS, Eq.(|34|), is its nonlin- 
ear dependence on the occupation number n, which leads 
to the switching effect and hysteresis in the I-V character- 
istics for d > 2, as will be shown below. It contains full 
information about all possible correlation and inelastic 
effects in transport, in particular, all the vibron-assisted 
tunneling processes (Fig. 2). 



IV. NONLINEAR RATE EQUATION AND 
SWITCHING 



a 




i — 











FIG. 2. Schematic of the inelastic processes assisting tun- 
neling through the molecular quantum dot (see Fig. 1). The 
tunneling through coupled electron-vibron system may pro- 
ceed with the emission (process e) or absorption (process a) 
of the vibrons. The absorption is possible only at non-zero 
temperatures. 



where ^1(2) are the transition rates from left (right) leads 
to MQD. Taking into account that p{uj) = d, Eq.@ 
for the symmetric leads, Ti = T2, reduces to 



2nd- J diup(Lu) (h + f 2 ). 



(36) 



which automatically satisfies < n < 1. Explicitly, the 
self-consistent equation for the occupation number is 



d-l 



= — Z r {n)[na r + (1 — n)b r 



(37) 



Generally, the electron density obeys an infinite set 
of rate equations for many-particle GFs which can be 
derived in the framework of a tunneling Hamiltonian in- 
cluding correlations jl6) . In the case of MQD only weakly 
coupled with leads one can apply the Fermi-Dirac golden 
rule to obtain an equation for n. Equating incoming and 
outgoing numbers of electrons in MQD per unit time we 
obtain the self-consistent equation for the level occupa- 
tion n as 

/DO 
du ! {T 1 f 1 {co) + r 2 f2^)}p^) 
-00 

/oo 
du {r x [1 - /i (w)] + T 2 [1 - / a (to)}} p(oj) = (35) 
-00 



where 



r = Z h (0 ( [fi(rU - lu ) + h{rU - Iuj q )} 



1=0 

+ (1 - 5 m )e-^[h{rU + luj ) + f 2 (rU + lu> )]\ , (38) 

OO / 

W = ZY, h (0 ( ^ [h{rU + lcu ) + f 2 (rU + Ilu )} 

+ (1 - S^e-^ihirU - ltu ) + f 2 (rU - lu> )}\ . (39) 
The current is expressed as 
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j = ^ = d J2 Z r(n)[na' r + (l-n)b' r ], 



(40) 



r=0 



P">o> 



where 

oo 

a' r = Z £ Ii (0 I [h(rU - lcu ) - f 2 (rU - lw )\ 
1=0 ^ 

+ (1 - d^e-^ 1 [h{rU + lco ) - f 2 {rU + Icj )} 

CO / 

b' r = Zj2ll (0 ( ^ [h(rU + M + h(rU + lu )] 
i=o ^ 

+ (1 - 8i )e-^ [h(rU - lw ) + f 2 (rU - lu )] 
Let us analyze the I-V curves for d = 1, 2, and 4. 



However, one of them is negative because < b r < a r < 1 
for any temperature and voltage. Therefore, we conclude 
that there is only one physical population of MQD, and 
the current is also single-valued at any voltage and tem- 
perature, in agreement with the Hubbard model [fl6|. 



(41) 



(42) 



A. Absence of switching for nondegenerate and 
two— fold degenerate MQD 

There is one term in the sum over r, r = with 
Zo(n) — 1, if c? = 1. Hence there is only one solution 
of the rate equation (37) 



bo 



2 + b - a 

and the current is single-valued at any voltage 



J 



26p + a' b - a Q b' 
2 + b - a 



(43) 



(44) 



This is an exact result, which is valid for any e-ph cou- 
pling and any phonon frequency. We have to conclude 
that there is no switching of a nondegenerate MQD. The 
opposite conclusion reached in Ref. |2^] might be due to 
the Born-Oppenheimer (static) approximation used by 
Gogolin and Komnik. In fact, the Born-Oppenheimer ap- 
proximation does not apply to the non-degenerate level 
model, since there are no "fast" , compared to phonon 
times l/wo, electron transitions within the "molecule". 

In the case of a double-degenerate MQD, d = 2, there 
are two terms, which contribute to the sum over r, with 
Zo(n) = l — rt and Z\(n) — n. The rate equation becomes 
a quadratic one 

n 2 (a + ax — b — &i) + n(2 - a + b a - b x ) — b = 0. 

(45) 

with two solutions, 

2 — a + b - bi 



ni.2 = 



± 



2(a + ax-bo - bx) 
(2-ao + bo-h) 2 



bo 



_4((a + ax - bo - bx) 2 a + ax - b - h 



1/2 



(46) 



B. Switching of four- fold degenerate MQD 



In this case the rate equation is of the fourth power in 



2n = (1 - n) 3 [na Q + (1 - n)b ] 

+3n(l - n) 2 [nax + (1 - n)bx] 
+3n 2 (l - n)[na 2 + (1 - n)b 2 ] 
+n 3 [na 3 + (1 - n)b 3 ]. 



(47) 



In the limit I7I 1 we have b r = a r , Z 
ing interaction is U = U c [sec Eq.([To| 
reduced to 

2n = (1 - nfb Q + 3n(l - n) 2 b x 
+3n 2 (l-n)b 2 +n 3 b 3 . 



■ 1, the remain- 
and Eq.(^|) is 

(48) 



If we assume now that the non-vibron interaction U c is 
negative, for example due to valence fluctuations, then 
we recover the negative-?/ model p6| , and the kinetic 
equation is reduced to 



2n = 1 - (1 - n) 3 



(49) 



in the voltage range A — \U\ < eV/2 < A and T = 
because bo = and b\ = b 2 — 63 = 1 there, if \U\ < A/2. 
The current is simplified as 



j = 2n. 



(50) 



Equation (|49|) has two physical roots, n = and n — 
(3 — 5 1 / 2 )/2 « 0.38. Hence we obtain two stationary 
states of MQD with low (zero at T = Q) and high cur- 
rent, j w 0.76 for the same voltage as we discussed earlier 
in Ref. [[l6). The current-voltage characteristics show a 
hysteretic behavior for d — 4. When the voltage increases 
from zero, 4-fold degenerate MQD remains in a low- 
current state until the threshold eV 2 /2 = A is reached. 
Remarkably, when the voltage decreases from the value 
above the threshold V 2 , the molecule remains in the high- 
current state down to the voltage eVi/2 = A — \U\ well 
below the threshold V 2 . This is a correlation mecha- 
nism of electronic molecular switching without retarda- 
tion. Therefore, the negative-?/ degenerate molecular dot 
possesses the volatile memory originating from the many- 
particle attractive correlations. 

The e-ph coupling results in the phonon sidebands of 
DOS, which are fully taken into account in Eq.(|35|) for the 
self-consistent occupation of the molecular level n. In the 
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FIG. 3. The I-V curves for tunneling through molecular 
quantum dot (Fig. 1) with the electron-vibron coupling con- 
stant ■y 2 = 11/13 and ujq/ A = 0.2. The up arrows show that 
the current picks up at some voltage when it is biased, and 
then drops at lower voltage when the bias is being reduced. 
The bias dependence of current basically repeats the shape of 
the level occupation n (right column). Steps on the curve cor- 
respond to the changing population of the phonon side-bands, 
which are shown in Fig. 1. The current hysteresis persist up 
to some critical temperature, which is low, T/A ~ 0.01. 



case of the Coulomb repulsion and the electron-vibron 
coupling the effective interaction will be attractive, if 



* 1 ^ 



(51) 



There is an important difference between switching 
with vibron-mediated electron-electron attraction and 
the negative-C/ model. We show the numerical results for 
uiq = 0.2 (in units of A, as all the energies in the prob- 
lem) and U c = for two values of the coupling constant, 
7 2 = 11/13 (Fig. 3) and 7 2 = 13/11 (Fig. 4). This 
case formally corresponds to U = —2 r y 2 uio «s —0.4, i.e. 
close to the same value of the attraction as we have used 
in the negative-l/ model (we selected those values 
of 7 2 to avoid accidental commensurability of a ladder 
of the correlated level energies separated by U and gen- 
erated by them phonon side-band ladders with the step 
uiq). In the negative-!/ case the threshold voltages were 
eV 1 /2A = 1 - 0.4 = 0.6 and eV 2 /2A = 1. However, in 
the vibron case the threshold for the onset of bistability 
corresponds to larger voltage bias compared to negative- 
U case (at eV/2A = 0.86 for 7 2 = 11/13 and lu = 0.2). 
The I-V in the vibron case is much more complex too: as 



FIG. 4. The I-V curves for tunneling through molecu- 
lar quantum dot with the electron-vibron coupling constant 
7 2 = 13/11, which is slightly larger than in Fig. 3. The I-V 
curves change substantially: the current pickup shifts to lower 
bias voltages and the curves show substantial change of shape 
compared to Fig. 3. The hysteresis persists to slightly higher 
temperatures, although also small, as compared to the previ- 
ous case. 



one can see from Figs. 3 and 4, the current discontinuity 
at the threshold strongly depends on the value of the e- 
ph coupling constant. The inelastic tunneling processes 
through the level, accompanied by emission/ absorption 
of the phonons (Fig. 2), manifest themselves as steps 
on the I-V curve, Figs. 3 and 4. Those steps are gen- 
erated by the phonon side-bands originating from corre- 
lated levels on the dot with the energies A, A + U, 
A + (d — 1)U. Since oj is not generally commensurate 
with U, we obtain pretty irregular picture of the steps on 
I-V curve. This comes as no surprise, since the kinetic 
equation ( [47| ) is much more complex compared to the one 
in the negative-U case, cf. (|9|). The level occupation ap- 
proaches the limiting value n = 0.5 at large bias voltages, 
Figs. 3 and 4. 

The bistability region shrinks down with temperature. 
In the specific example of the negative-C/ model with 
U/A = -0.4 the bistability is over at T/A w 0.1 [§J|. 
Importantly, in the vibron case this happens at much 
lower temperatures. Indeed, the hysteresis loop almost 
closes at T/A = 0.01. The critical temperature, below 
which the current bistability exists in the vibron case, is 
suppressed by about an order of magnitude compared to 
the negative-?/ case. At finite temperatures, the overall 
I-V curve shows the smoothed out steps at the bias volt- 
ages coinciding with the voltage where the vibrons are 
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emitted/ absorbed as at T = (Figs. 3, 4). 

In conclusion, we have developed the multi-polaron 
theory of tunneling through a molecular quantum dot 
(MQD) taking phonon sidebands and strong electron cor- 
relations into account. The degenerate MQD with strong 
electron- vibron coupling shows a hysteretic volatile mem- 
ory if the degeneracy of the molecular level is larger than 
two, d > 2. The hysteretic behavior strongly depends 
on electron- vibron coupling and characteristic vibron fre- 
quencies. The current bistability vanishes above some 
critical temperature. It would be very interesting to look 
for an experimental realization of the model, possibly in 
a system containing a certain conjugated central part, 
which exhibits the attractive correlations of carriers with 
large degeneracy d > 2. 
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